NOTES FROM SF:
FOR STEAMBOAT, COMPARE THE DETRENDED CLIAMTE NORMAL CURVES AGAINST EACH OTHER
Derry & fassnacht 2010 Compare climatology of snotel statuions= they don’t have the same climatol9ogy based on the snow data
Figure 5b - https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2009WR007835
This work is very preliminary as I get back into the coding swing of things. Data wrangling and figure generation will be done via R, but the rest of the project will be done using good ol’ microsoft products. This is just an entry point into data crunching and should by no means be considered a final product.
Also, I’m not great at this but whatever. I could automate this, but I’ll figure that out shortly!
Need to update figures & analyze sites for seasonal trends
SNOTEL data was provided by the NRCS. Data was cleaned by removing outliers that are likely implausible; any year with more than 15 observations missing was removed. Temperatures were adjusted using the Morrisey method for stations identified by Ma et al (2019) due to SNOTEL temperature sensor changes, with the adjustment applied to pre-sensor change data. Daily mean observations were detrended to determine whether values were increasing or decreasing from the entire time series trend. Daily mean temperatures were first averaged by water year, with all water year means then averaged by day of water year. The mean temperature by day for the period of record was averaged. To find the standard deviation, the daily mean temperatures by water year was subtracted from the averaged mean temperature by day for the period of record. All water year means averaged by day of water year were subtracted from the temperature mean. The resulting values were then added together to find the “residual” of the daily mean temperatures by water year. The standard deviation was then computed from those residuals, with trends analyzed by Mann‐Kendall significance test and Theil‐Sen’s rate of change. Significant trends are identified with p-values of less than 0.10.
Morrisey Method
The Morrisey Method is taken from Ma, Fassnacht and Kampf..
In R script: T(adjusted) = 5.3x10(-7)xT(old)4+3.72x10(-5)xT(old)3-2.16x10(-3)xT(old)2-7.32x10^(-2)xT(old)+1.37
In the Ma et al. spreadsheet, H1 is Morrisey, H2 is Oiler
Yampa Area SNOTEL sites:
Burro Mountain 378 NSCE- Original, Bias- Morrisey 7/11/2002 (NSCE- 0.87 vs 0.84, Bias- -0.03 vs. 0.11)
Columbine 408 Morrisey 7/22/2005
Columbine Pass 409 Morrisey 6/23/2005
Crosho 426 Morrisey 7/21/2005
Dry Lake 457 Oyler -> Morrisey 7/30/2003
Elk River 467 Oyler -> Morrisey 8/7/2006
Lynx Pass 607 Morrisey 5/22/2006
Rabbit Ears 709 Morrisey 8/7/2006 (Does not go.)
Ripple Creek 717 Oyler -> Morrisey 8/7/2006
Tower 825 Original 8/18/2004
Trapper Lake 827 Original 12/13/2004
SNOTEL_yampa_area <- snotel_download(site_id = c(378, 408, 409, 426, 457, 467, 607, 709, 717, 825, 827), path = tempdir('../data'), internal = TRUE)
write.csv(SNOTEL_yampa_area,"C:/Users/13074/Documents/ESS580/thesis_project/Yampa/data_raw/snotel_yampa.csv", row.names = FALSE) #write in the raw data
SNOTEL_yampa_area <- read.csv("C:/Users/13074/Documents/ESS580/thesis_project/Yampa/data_raw/snotel_yampa.csv", header = TRUE)
NSCE- Original, Bias- Morrisey 7/11/2002 **Original is marked as having the smallest bias, but the sheet says original is smaller.)
snotel_378 <- SNOTEL_yampa_area %>%
filter(site_id == "378")
#str(snotel_378) # check the date, usually a character.
snotel_378$Date <- as.Date(snotel_378$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.
#THIS WILL CHANGE FOR EACH STATION
snotel_378_clean <- snotel_378 %>% # filter for the timeframe
filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
#filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers
addWaterYear() %>%
mutate(daymonth = format(as.Date(Date), "%d-%m")) %>%
na.omit()
#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))
snotel_378_clean <- snotel_378_clean %>%
group_by(waterYear)%>%
mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))
# Check for outliers
ggplot(snotel_378_clean, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
snotel_378_clean <- snotel_378_clean %>%
mutate(temp_diff = abs(temperature_min - temperature_max)) %>%
filter(temperature_mean > -40)# %>%
#filter(temperature_mean < 25)
ggplot(snotel_378_clean, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”
# filtering for temperature anomalies
#snotel_378_cull_count <- snotel_378_clean %>%
# filter(temperature_min > -40) %>%
# count(waterYear)
#snotel_378_cull_count
# filtering for too few observations in a year
snotel_378_cull_count_days <- snotel_378_clean %>%
group_by(waterYear) %>%
count(waterYear) %>%
filter(n < 350)
snotel_378_cull_count_days
## # A tibble: 10 x 2
## # Groups: waterYear [10]
## waterYear n
## <dbl> <int>
## 1 1985 1
## 2 1988 298
## 3 1989 341
## 4 1993 281
## 5 1994 261
## 6 1995 344
## 7 1996 320
## 8 2002 312
## 9 2016 310
## 10 2022 338
snotel_378_clean_culled <- snotel_378_clean %>%
filter(waterYear != "1985" & waterYear != "1986" & waterYear != "1988" & waterYear != "1989" & waterYear != "1993" & waterYear != "1994" & waterYear != "1995" & waterYear != "1996" & waterYear != "2002" & waterYear != "2016" & waterYear != "2022")# & waterYear != "2017") #%>%
#filter(temperature_mean > -49)
ggplot(snotel_378_clean_culled, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
Culling WY 1986 as well.
ggplot(snotel_378_clean_culled, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
temp_378_xts <- xts(snotel_378_clean_culled$temperature_mean, order.by = snotel_378_clean_culled$Date)
dygraph(temp_378_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)")
#snotel_378_clean_culled <- snotel_378_clean_culled %>%
# filter(temperature_mean < 19.3)
temp_378_xts <- xts(snotel_378_clean_culled$temperature_mean, order.by = snotel_378_clean_culled$Date)
dygraph(temp_378_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)")
NSCE- Original, Bias- Morrisey 7/11/2002
snotel_378_adjusted <- snotel_378_clean_culled %>%
mutate(temp_ad = if_else(Date < "2002-07-11", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))
Non-corrected
#using the clean culled df:
#average water year temperature
yearly_wy_aver_378 <- snotel_378_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp = mean(temperature_mean))
#Average temperature by day for all water years:
daily_wy_aver_378 <- yearly_wy_aver_378 %>%
group_by(daymonth) %>%
mutate(aver_day_temp = mean(aver_ann_temp))
#average mean temperature by day for the period of record:
daily_wy_aver_378 <- daily_wy_aver_378 %>%
group_by(daymonth) %>%
mutate(all_ave_temp = mean(daily_wy_aver_378$aver_day_temp))
# try to show all years as means.
daily_wy_aver2_378 <-daily_wy_aver_378 %>%
group_by(waterDay) %>%
mutate(date_temp = mean(temperature_mean))
daily_wy_aver2_378$date_temp <- signif(daily_wy_aver2_378$date_temp,3) #reduce the sig figs
ggplot(daily_wy_aver2_378, aes(x = waterDay, y = date_temp))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
standard_dev_378 <- daily_wy_aver_378 %>%
group_by(waterYear) %>%
mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_378 <- standard_dev_378 %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_378 <- standard_dev_all_378 %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_378 %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')
| waterYear | sd_2 |
|---|---|
| 1987 | 8.562568 |
| 1990 | 9.114355 |
| 1991 | 9.254156 |
| 1992 | 8.275786 |
| 1997 | 8.898932 |
| 1998 | 8.983050 |
| 1999 | 8.073327 |
| 2000 | 8.686540 |
| 2001 | 9.301471 |
| 2003 | 8.509538 |
| 2004 | 8.017044 |
| 2005 | 7.723002 |
| 2006 | 8.547829 |
| 2007 | 8.796189 |
| 2008 | 8.749038 |
| 2009 | 7.771973 |
| 2010 | 8.617066 |
| 2011 | 8.282002 |
| 2012 | 8.445614 |
| 2013 | 8.991985 |
| 2014 | 8.166025 |
| 2015 | 7.657629 |
| 2017 | 8.034586 |
| 2018 | 8.020621 |
| 2019 | 8.518154 |
| 2020 | 8.848254 |
| 2021 | 8.718771 |
ggplot(standard_dev_all_378, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Non-corrected standard deviation of SNOTEL 378 average temperatures for water years 2005-2021
sd_mk_378 <- mk.test(standard_dev_all_378$sd_2)
print(sd_mk_378)
##
## Mann-Kendall trend test
##
## data: standard_dev_all_378$sd_2
## z = -1.5427, n = 27, p-value = 0.1229
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## -75.0000000 2301.0000000 -0.2136752
sd_sens_378 <- sens.slope(standard_dev_all_378$sd_2)
print(sd_sens_378)
##
## Sen's slope
##
## data: standard_dev_all_378$sd_2
## z = -1.5427, n = 27, p-value = 0.1229
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.047456119 0.006621254
## sample estimates:
## Sen's slope
## -0.0230883
#using the clean culled df:
#average water year temperature
yearly_wy_aver_378_ad <- snotel_378_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_378_ad <- yearly_wy_aver_378_ad %>%
group_by(daymonth) %>%
mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_378_ad <- daily_wy_aver_378_ad %>%
group_by(daymonth) %>%
mutate(all_ave_temp_ad = mean(daily_wy_aver_378_ad$aver_day_temp_ad))
# try to show all years as means.
daily_wy_aver2_378_ad <-daily_wy_aver_378_ad %>%
group_by(waterDay) %>%
mutate(date_temp_ad = mean(temp_ad))
daily_wy_aver2_378_ad$date_temp_ad <- signif(daily_wy_aver2_378_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_378_ad, aes(x = waterDay, y = date_temp_ad))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
standard_dev_378_ad <- daily_wy_aver_378_ad %>%
group_by(waterYear) %>%
#filter(waterYear >= 1987 & waterYear <= 2021) %>%
mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_378_ad <- standard_dev_378_ad %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_378_ad <- standard_dev_all_378_ad %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_378_ad %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')
| waterYear | sd_2 |
|---|---|
| 1987 | 7.940547 |
| 1990 | 8.436056 |
| 1991 | 8.633602 |
| 1992 | 7.670874 |
| 1997 | 8.260208 |
| 1998 | 8.299486 |
| 1999 | 7.471734 |
| 2000 | 7.993287 |
| 2001 | 8.584377 |
| 2003 | 8.508522 |
| 2004 | 8.016141 |
| 2005 | 7.722598 |
| 2006 | 8.547499 |
| 2007 | 8.795362 |
| 2008 | 8.748133 |
| 2009 | 7.771479 |
| 2010 | 8.616338 |
| 2011 | 8.281602 |
| 2012 | 8.444783 |
| 2013 | 8.991116 |
| 2014 | 8.165236 |
| 2015 | 7.657576 |
| 2017 | 8.034186 |
| 2018 | 8.019986 |
| 2019 | 8.517552 |
| 2020 | 8.847387 |
| 2021 | 8.717958 |
ggplot(standard_dev_all_378_ad, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Morrisey-corrected standard deviation of SNOTEL 378 average temperatures for water years 1986-2021
sd_mk_378_ad <- mk.test(standard_dev_all_378_ad$sd_2)
print(sd_mk_378_ad)
##
## Mann-Kendall trend test
##
## data: standard_dev_all_378_ad$sd_2
## z = 1.3342, n = 27, p-value = 0.1821
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## 65.0000000 2301.0000000 0.1851852
sd_sens_378_ad <- sens.slope(standard_dev_all_378_ad$sd_2)
print(sd_sens_378_ad)
##
## Sen's slope
##
## data: standard_dev_all_378_ad$sd_2
## z = 1.3342, n = 27, p-value = 0.1821
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.009653315 0.035911433
## sample estimates:
## Sen's slope
## 0.01222034
Morrisey 7/22/2005
snotel_408 <- SNOTEL_yampa_area %>%
filter(site_id == "408")
#str(snotel_408) # check the date, usually a character.
snotel_408$Date <- as.Date(snotel_408$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.
#THIS WILL CHANGE FOR EACH STATION
snotel_408_clean <- snotel_408 %>% # filter for the timeframe
filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
#filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers
addWaterYear() %>%
mutate(daymonth = format(as.Date(Date), "%d-%m")) %>%
na.omit()
#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))
snotel_408_clean <- snotel_408_clean %>%
group_by(waterYear)%>%
mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))
# Check for outliers
ggplot(snotel_408_clean, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
snotel_408_clean <- snotel_408_clean %>%
mutate(temp_diff = abs(temperature_min - temperature_max)) %>%
filter(temperature_mean > -40)# %>%
#filter(temperature_mean < 25)
ggplot(snotel_408_clean, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”
# filtering for temperature anomalies
#snotel_408_cull_count <- snotel_408_clean %>%
# filter(temperature_min > -40) %>%
# count(waterYear)
#snotel_408_cull_count
# filtering for too few observations in a year
snotel_408_cull_count_days <- snotel_408_clean %>%
group_by(waterYear) %>%
count(waterYear) %>%
filter(n < 350)
snotel_408_cull_count_days
## # A tibble: 8 x 2
## # Groups: waterYear [8]
## waterYear n
## <dbl> <int>
## 1 1981 343
## 2 1983 317
## 3 1984 234
## 4 1985 195
## 5 1986 336
## 6 1987 298
## 7 1994 272
## 8 2002 347
snotel_408_clean_culled <- snotel_408_clean %>%
filter(waterYear > "1987" & waterYear != "1993" & waterYear != "1994" & waterYear != "2002")# & waterYear != "1985" & waterYear != "1986" & waterYear != "1987" & waterYear != "1994" & waterYear != "2002")# & waterYear != "2002" & waterYear != "2016" & waterYear != "2022")# & waterYear != "2017") #%>%
#filter(temperature_mean > -49)
ggplot(snotel_408_clean_culled, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
Culling WY 1986 & 1993 as well.
ggplot(snotel_408_clean_culled, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
temp_408_xts <- xts(snotel_408_clean_culled$temperature_mean, order.by = snotel_408_clean_culled$Date)
dygraph(temp_408_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)")
snotel_408_clean_culled <- snotel_408_clean_culled %>%
filter(temperature_mean > -30)
temp_408_xts <- xts(snotel_408_clean_culled$temperature_mean, order.by = snotel_408_clean_culled$Date)
dygraph(temp_408_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)")
snotel_408_adjusted <- snotel_408_clean_culled %>%
mutate(temp_ad = if_else(Date < "2005-07-22", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))
Non-corrected
#using the clean culled df:
#average water year temperature
yearly_wy_aver_408 <- snotel_408_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp = mean(temperature_mean))
#Average temperature by day for all water years:
daily_wy_aver_408 <- yearly_wy_aver_408 %>%
group_by(daymonth) %>%
mutate(aver_day_temp = mean(aver_ann_temp))
#average mean temperature by day for the period of record:
daily_wy_aver_408 <- daily_wy_aver_408 %>%
group_by(daymonth) %>%
mutate(all_ave_temp = mean(daily_wy_aver_408$aver_day_temp))
# try to show all years as means.
daily_wy_aver2_408 <-daily_wy_aver_408 %>%
group_by(waterDay) %>%
mutate(date_temp = mean(temperature_mean))
daily_wy_aver2_408$date_temp <- signif(daily_wy_aver2_408$date_temp,3) #reduce the sig figs
ggplot(daily_wy_aver2_408, aes(x = waterDay, y = date_temp))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
standard_dev_408 <- daily_wy_aver_408 %>%
group_by(waterYear) %>%
mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_408 <- standard_dev_408 %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_408 <- standard_dev_all_408 %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_408 %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')
| waterYear | sd_2 |
|---|---|
| 1988 | 9.911311 |
| 1989 | 9.479817 |
| 1990 | 8.845588 |
| 1991 | 9.141713 |
| 1992 | 8.235717 |
| 1995 | 8.387811 |
| 1996 | 8.807923 |
| 1997 | 8.867338 |
| 1998 | 9.074387 |
| 1999 | 8.232888 |
| 2000 | 8.585091 |
| 2001 | 9.380546 |
| 2003 | 8.847298 |
| 2004 | 8.495495 |
| 2005 | 8.135796 |
| 2006 | 8.734304 |
| 2007 | 8.839901 |
| 2008 | 8.888992 |
| 2009 | 7.926154 |
| 2010 | 8.738675 |
| 2011 | 8.455577 |
| 2012 | 8.702054 |
| 2013 | 9.110505 |
| 2014 | 8.268271 |
| 2015 | 7.870630 |
| 2016 | 8.492743 |
| 2017 | 8.238399 |
| 2018 | 8.140089 |
| 2019 | 8.644103 |
| 2020 | 9.027143 |
| 2021 | 8.902353 |
| 2022 | 8.706784 |
ggplot(standard_dev_all_408, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Non-corrected standard deviation of SNOTEL 408 average temperatures for water years 2005-2021
sd_mk_408 <- mk.test(standard_dev_all_408$sd_2)
print(sd_mk_408)
##
## Mann-Kendall trend test
##
## data: standard_dev_all_408$sd_2
## z = -1.7676, n = 32, p-value = 0.07713
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## -110.0000000 3802.6666667 -0.2217742
sd_sens_408 <- sens.slope(standard_dev_all_408$sd_2)
print(sd_sens_408)
##
## Sen's slope
##
## data: standard_dev_all_408$sd_2
## z = -1.7676, n = 32, p-value = 0.07713
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.033742130 0.001027723
## sample estimates:
## Sen's slope
## -0.01615257
#using the clean culled df:
#average water year temperature
yearly_wy_aver_408_ad <- snotel_408_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_408_ad <- yearly_wy_aver_408_ad %>%
group_by(daymonth) %>%
mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_408_ad <- daily_wy_aver_408_ad %>%
group_by(daymonth) %>%
mutate(all_ave_temp_ad = mean(daily_wy_aver_408_ad$aver_day_temp_ad))
# try to show all years as means.
daily_wy_aver2_408_ad <-daily_wy_aver_408_ad %>%
group_by(waterDay) %>%
mutate(date_temp_ad = mean(temp_ad))
daily_wy_aver2_408_ad$date_temp_ad <- signif(daily_wy_aver2_408_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_408_ad, aes(x = waterDay, y = date_temp_ad))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
standard_dev_408_ad <- daily_wy_aver_408_ad %>%
group_by(waterYear) %>%
#filter(waterYear >= 1987 & waterYear <= 2021) %>%
mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_408_ad <- standard_dev_408_ad %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_408_ad <- standard_dev_all_408_ad %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_408_ad %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')
| waterYear | sd_2 |
|---|---|
| 1988 | 9.255004 |
| 1989 | 8.873391 |
| 1990 | 8.216764 |
| 1991 | 8.545580 |
| 1992 | 7.664569 |
| 1995 | 7.780528 |
| 1996 | 8.186934 |
| 1997 | 8.260524 |
| 1998 | 8.410419 |
| 1999 | 7.651225 |
| 2000 | 7.928801 |
| 2001 | 8.691699 |
| 2003 | 8.174582 |
| 2004 | 7.895518 |
| 2005 | 7.482674 |
| 2006 | 8.735962 |
| 2007 | 8.841502 |
| 2008 | 8.890875 |
| 2009 | 7.928502 |
| 2010 | 8.741101 |
| 2011 | 8.457260 |
| 2012 | 8.703940 |
| 2013 | 9.112301 |
| 2014 | 8.270874 |
| 2015 | 7.872234 |
| 2016 | 8.495258 |
| 2017 | 8.239886 |
| 2018 | 8.142322 |
| 2019 | 8.646452 |
| 2020 | 9.028953 |
| 2021 | 8.904182 |
| 2022 | 8.708980 |
ggplot(standard_dev_all_408_ad, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Morrisey-corrected standard deviation of SNOTEL 408 average temperatures for water years 1986-2021
sd_mk_408_ad <- mk.test(standard_dev_all_408_ad$sd_2)
print(sd_mk_408_ad)
##
## Mann-Kendall trend test
##
## data: standard_dev_all_408_ad$sd_2
## z = 1.2487, n = 32, p-value = 0.2118
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## 78.0000000 3802.6666667 0.1572581
sd_sens_408_ad <- sens.slope(standard_dev_all_408_ad$sd_2)
print(sd_sens_408_ad)
##
## Sen's slope
##
## data: standard_dev_all_408_ad$sd_2
## z = 1.2487, n = 32, p-value = 0.2118
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.007349474 0.033195852
## sample estimates:
## Sen's slope
## 0.01316118
Morrisey 6/23/2005
snotel_409 <- SNOTEL_yampa_area %>%
filter(site_id == "409")
#str(snotel_409) # check the date, usually a character.
snotel_409$Date <- as.Date(snotel_409$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.
#THIS WILL CHANGE FOR EACH STATION
snotel_409_clean <- snotel_409 %>% # filter for the timeframe
filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
#filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers
addWaterYear() %>%
mutate(daymonth = format(as.Date(Date), "%d-%m")) %>%
na.omit()
#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))
snotel_409_clean <- snotel_409_clean %>%
group_by(waterYear)%>%
mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))
# Check for outliers
ggplot(snotel_409_clean, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
snotel_409_clean <- snotel_409_clean %>%
mutate(temp_diff = abs(temperature_min - temperature_max)) %>%
filter(temperature_mean > -50) %>%
filter(temperature_mean < 40)
ggplot(snotel_409_clean, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”
# filtering for temperature anomalies
#snotel_409_cull_count <- snotel_409_clean %>%
# filter(temperature_min > -40) %>%
# count(waterYear)
#snotel_409_cull_count
# filtering for too few observations in a year
snotel_409_cull_count_days <- snotel_409_clean %>%
group_by(waterYear) %>%
count(waterYear) %>%
filter(n < 350)
snotel_409_cull_count_days
## # A tibble: 5 x 2
## # Groups: waterYear [5]
## waterYear n
## <dbl> <int>
## 1 1994 331
## 2 1995 190
## 3 2005 329
## 4 2019 348
## 5 2022 346
snotel_409_clean_culled <- snotel_409_clean %>%
filter(waterYear != "1994" & waterYear != "1995" & waterYear != "2005" & waterYear != "2019" & waterYear != "2022")# & waterYear != "1986" & waterYear != "1987" & waterYear != "1994" & waterYear != "2002")# & waterYear != "2002" & waterYear != "2016" & waterYear != "2022")# & waterYear != "2017") #%>%
#filter(temperature_mean > -49)
ggplot(snotel_409_clean_culled, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
ggplot(snotel_409_clean_culled, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
temp_409_xts <- xts(snotel_409_clean_culled$temperature_mean, order.by = snotel_409_clean_culled$Date)
dygraph(temp_409_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)")
#snotel_409_clean_culled <- snotel_409_clean_culled %>%
# filter(temperature_mean > -30)
#temp_409_xts <- xts(snotel_409_clean_culled$temperature_mean, order.by = snotel_409_clean_culled$Date)
#dygraph(temp_409_xts) %>%
# dyAxis("y", label = "Daily mean temperature (°C)")
Morrisey 6/23/2005
snotel_409_adjusted <- snotel_409_clean_culled %>%
mutate(temp_ad = if_else(Date < "2005-06-23", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))
Non-corrected
#using the clean culled df:
#average water year temperature
yearly_wy_aver_409 <- snotel_409_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp = mean(temperature_mean))
#Average temperature by day for all water years:
daily_wy_aver_409 <- yearly_wy_aver_409 %>%
group_by(daymonth) %>%
mutate(aver_day_temp = mean(aver_ann_temp))
#average mean temperature by day for the period of record:
daily_wy_aver_409 <- daily_wy_aver_409 %>%
group_by(daymonth) %>%
mutate(all_ave_temp = mean(daily_wy_aver_409$aver_day_temp))
# try to show all years as means.
daily_wy_aver2_409 <-daily_wy_aver_409 %>%
group_by(waterDay) %>%
mutate(date_temp = mean(temperature_mean))
daily_wy_aver2_409$date_temp <- signif(daily_wy_aver2_409$date_temp,3) #reduce the sig figs
ggplot(daily_wy_aver2_409, aes(x = waterDay, y = date_temp))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
standard_dev_409 <- daily_wy_aver_409 %>%
group_by(waterYear) %>%
mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_409 <- standard_dev_409 %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_409 <- standard_dev_all_409 %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_409 %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')
| waterYear | sd_2 |
|---|---|
| 1987 | 8.326403 |
| 1988 | 9.167257 |
| 1989 | 8.983106 |
| 1990 | 8.555010 |
| 1991 | 7.839925 |
| 1992 | 7.987351 |
| 1993 | 8.471182 |
| 1996 | 8.508682 |
| 1997 | 8.478813 |
| 1998 | 8.929006 |
| 1999 | 7.623863 |
| 2000 | 8.631812 |
| 2001 | 9.047367 |
| 2002 | 9.516272 |
| 2003 | 8.791318 |
| 2004 | 8.708583 |
| 2006 | 8.125432 |
| 2007 | 8.536611 |
| 2008 | 8.816576 |
| 2009 | 7.846243 |
| 2010 | 8.511418 |
| 2011 | 8.496243 |
| 2012 | 8.285134 |
| 2013 | 8.792304 |
| 2014 | 7.958185 |
| 2015 | 7.393468 |
| 2016 | 8.187794 |
| 2017 | 8.066498 |
| 2018 | 7.805356 |
| 2020 | 8.856536 |
| 2021 | 8.674169 |
ggplot(standard_dev_all_409, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Non-corrected standard deviation of SNOTEL 409 average temperatures for water years 2005-2021
sd_mk_409 <- mk.test(standard_dev_all_409$sd_2)
print(sd_mk_409)
##
## Mann-Kendall trend test
##
## data: standard_dev_all_409$sd_2
## z = -1.1218, n = 31, p-value = 0.262
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## -67.000000 3461.666667 -0.144086
sd_sens_409 <- sens.slope(standard_dev_all_409$sd_2)
print(sd_sens_409)
##
## Sen's slope
##
## data: standard_dev_all_409$sd_2
## z = -1.1218, n = 31, p-value = 0.262
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.034898600 0.009250748
## sample estimates:
## Sen's slope
## -0.01256069
#using the clean culled df:
#average water year temperature
yearly_wy_aver_409_ad <- snotel_409_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_409_ad <- yearly_wy_aver_409_ad %>%
group_by(daymonth) %>%
mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_409_ad <- daily_wy_aver_409_ad %>%
group_by(daymonth) %>%
mutate(all_ave_temp_ad = mean(daily_wy_aver_409_ad$aver_day_temp_ad))
# try to show all years as means.
daily_wy_aver2_409_ad <-daily_wy_aver_409_ad %>%
group_by(waterDay) %>%
mutate(date_temp_ad = mean(temp_ad))
daily_wy_aver2_409_ad$date_temp_ad <- signif(daily_wy_aver2_409_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_409_ad, aes(x = waterDay, y = date_temp_ad))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
standard_dev_409_ad <- daily_wy_aver_409_ad %>%
group_by(waterYear) %>%
#filter(waterYear >= 1987 & waterYear <= 2021) %>%
mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_409_ad <- standard_dev_409_ad %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_409_ad <- standard_dev_all_409_ad %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_409_ad %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')
| waterYear | sd_2 |
|---|---|
| 1987 | 7.660579 |
| 1988 | 8.456632 |
| 1989 | 8.306699 |
| 1990 | 7.862004 |
| 1991 | 7.209890 |
| 1992 | 7.340956 |
| 1993 | 7.806529 |
| 1996 | 7.809598 |
| 1997 | 7.810012 |
| 1998 | 8.200421 |
| 1999 | 7.004624 |
| 2000 | 7.910753 |
| 2001 | 8.302261 |
| 2002 | 8.749019 |
| 2003 | 8.046898 |
| 2004 | 8.036896 |
| 2006 | 8.125530 |
| 2007 | 8.536428 |
| 2008 | 8.816453 |
| 2009 | 7.846083 |
| 2010 | 8.511060 |
| 2011 | 8.495994 |
| 2012 | 8.285156 |
| 2013 | 8.792015 |
| 2014 | 7.958138 |
| 2015 | 7.393110 |
| 2016 | 8.187413 |
| 2017 | 8.066215 |
| 2018 | 7.805282 |
| 2020 | 8.856204 |
| 2021 | 8.673923 |
ggplot(standard_dev_all_409_ad, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Morrisey-corrected standard deviation of SNOTEL 409 average temperatures for water years 1986-2021
sd_mk_409_ad <- mk.test(standard_dev_all_409_ad$sd_2)
print(sd_mk_409_ad)
##
## Mann-Kendall trend test
##
## data: standard_dev_all_409_ad$sd_2
## z = 2.1415, n = 31, p-value = 0.03223
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## 127.0000000 3461.6666667 0.2731183
sd_sens_409_ad <- sens.slope(standard_dev_all_409_ad$sd_2)
print(sd_sens_409_ad)
##
## Sen's slope
##
## data: standard_dev_all_409_ad$sd_2
## z = 2.1415, n = 31, p-value = 0.03223
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## 0.002607708 0.039439754
## sample estimates:
## Sen's slope
## 0.02152587
Morrisey 7/21/2005
snotel_426 <- SNOTEL_yampa_area %>%
filter(site_id == "426")
#str(snotel_426) # check the date, usually a character.
snotel_426$Date <- as.Date(snotel_426$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.
#THIS WILL CHANGE FOR EACH STATION
snotel_426_clean <- snotel_426 %>% # filter for the timeframe
filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
#filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers
addWaterYear() %>%
mutate(daymonth = format(as.Date(Date), "%d-%m")) %>%
na.omit()
#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))
snotel_426_clean <- snotel_426_clean %>%
group_by(waterYear)%>%
mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))
# Check for outliers
ggplot(snotel_426_clean, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
snotel_426_clean <- snotel_426_clean %>%
mutate(temp_diff = abs(temperature_min - temperature_max)) %>%
filter(temperature_mean > -50) %>%
filter(temperature_mean < 40)
ggplot(snotel_426_clean, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”
# filtering for temperature anomalies
#snotel_426_cull_count <- snotel_426_clean %>%
# filter(temperature_min > -40) %>%
# count(waterYear)
#snotel_426_cull_count
# filtering for too few observations in a year
snotel_426_cull_count_days <- snotel_426_clean %>%
group_by(waterYear) %>%
count(waterYear) %>%
filter(n < 350)
snotel_426_cull_count_days
## # A tibble: 2 x 2
## # Groups: waterYear [2]
## waterYear n
## <dbl> <int>
## 1 2009 340
## 2 2021 349
snotel_426_clean_culled <- snotel_426_clean %>%
filter(waterYear != "2009" & waterYear != "2021")# & waterYear != "2005" & waterYear != "2019" & waterYear != "2022")# & waterYear != "1986" & waterYear != "1987" & waterYear != "1994" & waterYear != "2002")# & waterYear != "2002" & waterYear != "2016" & waterYear != "2022")# & waterYear != "2017") #%>%
#filter(temperature_mean > -49)
ggplot(snotel_426_clean_culled, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
ggplot(snotel_426_clean_culled, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
temp_426_xts <- xts(snotel_426_clean_culled$temperature_mean, order.by = snotel_426_clean_culled$Date)
dygraph(temp_426_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)")
#snotel_426_clean_culled <- snotel_426_clean_culled %>%
# filter(temperature_mean > -30)
#temp_426_xts <- xts(snotel_426_clean_culled$temperature_mean, order.by = snotel_426_clean_culled$Date)
#dygraph(temp_426_xts) %>%
# dyAxis("y", label = "Daily mean temperature (°C)")
snotel_426_adjusted <- snotel_426_clean_culled %>%
mutate(temp_ad = if_else(Date < "2005-07-21", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))
Non-corrected
#using the clean culled df:
#average water year temperature
yearly_wy_aver_426 <- snotel_426_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp = mean(temperature_mean))
#Average temperature by day for all water years:
daily_wy_aver_426 <- yearly_wy_aver_426 %>%
group_by(daymonth) %>%
mutate(aver_day_temp = mean(aver_ann_temp))
#average mean temperature by day for the period of record:
daily_wy_aver_426 <- daily_wy_aver_426 %>%
group_by(daymonth) %>%
mutate(all_ave_temp = mean(daily_wy_aver_426$aver_day_temp))
# try to show all years as means.
daily_wy_aver2_426 <-daily_wy_aver_426 %>%
group_by(waterDay) %>%
mutate(date_temp = mean(temperature_mean))
daily_wy_aver2_426$date_temp <- signif(daily_wy_aver2_426$date_temp,3) #reduce the sig figs
ggplot(daily_wy_aver2_426, aes(x = waterDay, y = date_temp))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
standard_dev_426 <- daily_wy_aver_426 %>%
group_by(waterYear) %>%
mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_426 <- standard_dev_426 %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_426 <- standard_dev_all_426 %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_426 %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')
| waterYear | sd_2 |
|---|---|
| 1987 | 8.642717 |
| 1988 | 8.008775 |
| 1989 | 9.449036 |
| 1990 | 8.683053 |
| 1991 | 9.170577 |
| 1992 | 8.537791 |
| 1993 | 8.676514 |
| 1994 | 9.525407 |
| 1995 | 8.192424 |
| 1996 | 8.622041 |
| 1997 | 8.670823 |
| 1998 | 8.768147 |
| 1999 | 7.907919 |
| 2000 | 8.392275 |
| 2001 | 9.164446 |
| 2002 | 9.813834 |
| 2003 | 8.961799 |
| 2004 | 8.571169 |
| 2005 | 8.096051 |
| 2006 | 8.582421 |
| 2007 | 8.845773 |
| 2008 | 8.787455 |
| 2010 | 8.872551 |
| 2011 | 8.348895 |
| 2012 | 8.502421 |
| 2013 | 8.997153 |
| 2014 | 8.050838 |
| 2015 | 7.657484 |
| 2016 | 8.378217 |
| 2017 | 8.126881 |
| 2018 | 8.079436 |
| 2019 | 8.651915 |
| 2020 | 8.964956 |
| 2022 | 8.545430 |
ggplot(standard_dev_all_426, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Non-corrected standard deviation of SNOTEL 426 average temperatures for water years 2005-2021
sd_mk_426 <- mk.test(standard_dev_all_426$sd_2)
print(sd_mk_426)
##
## Mann-Kendall trend test
##
## data: standard_dev_all_426$sd_2
## z = -1.2749, n = 34, p-value = 0.2023
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## -87.0000000 4550.3333333 -0.1550802
sd_sens_426 <- sens.slope(standard_dev_all_426$sd_2)
print(sd_sens_426)
##
## Sen's slope
##
## data: standard_dev_all_426$sd_2
## z = -1.2749, n = 34, p-value = 0.2023
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.02882424 0.00714949
## sample estimates:
## Sen's slope
## -0.009821212
#using the clean culled df:
#average water year temperature
yearly_wy_aver_426_ad <- snotel_426_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_426_ad <- yearly_wy_aver_426_ad %>%
group_by(daymonth) %>%
mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_426_ad <- daily_wy_aver_426_ad %>%
group_by(daymonth) %>%
mutate(all_ave_temp_ad = mean(daily_wy_aver_426_ad$aver_day_temp_ad))
# try to show all years as means.
daily_wy_aver2_426_ad <-daily_wy_aver_426_ad %>%
group_by(waterDay) %>%
mutate(date_temp_ad = mean(temp_ad))
daily_wy_aver2_426_ad$date_temp_ad <- signif(daily_wy_aver2_426_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_426_ad, aes(x = waterDay, y = date_temp_ad))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
standard_dev_426_ad <- daily_wy_aver_426_ad %>%
group_by(waterYear) %>%
#filter(waterYear >= 1987 & waterYear <= 2021) %>%
mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_426_ad <- standard_dev_426_ad %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_426_ad <- standard_dev_all_426_ad %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_426_ad %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')
| waterYear | sd_2 |
|---|---|
| 1987 | 8.032449 |
| 1988 | 7.469399 |
| 1989 | 8.830085 |
| 1990 | 8.043543 |
| 1991 | 8.571182 |
| 1992 | 7.909178 |
| 1993 | 8.068779 |
| 1994 | 8.814518 |
| 1995 | 7.580772 |
| 1996 | 7.973712 |
| 1997 | 8.044413 |
| 1998 | 8.100725 |
| 1999 | 7.332358 |
| 2000 | 7.718765 |
| 2001 | 8.465805 |
| 2002 | 9.079917 |
| 2003 | 8.257697 |
| 2004 | 7.940383 |
| 2005 | 7.444764 |
| 2006 | 8.582735 |
| 2007 | 8.846403 |
| 2008 | 8.788093 |
| 2010 | 8.872827 |
| 2011 | 8.349023 |
| 2012 | 8.502544 |
| 2013 | 8.997369 |
| 2014 | 8.051074 |
| 2015 | 7.658079 |
| 2016 | 8.378828 |
| 2017 | 8.126888 |
| 2018 | 8.080106 |
| 2019 | 8.652534 |
| 2020 | 8.965133 |
| 2022 | 8.546170 |
ggplot(standard_dev_all_426_ad, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Morrisey-corrected standard deviation of SNOTEL 426 average temperatures for water years 1986-2021
sd_mk_426_ad <- mk.test(standard_dev_all_426_ad$sd_2)
print(sd_mk_426_ad)
##
## Mann-Kendall trend test
##
## data: standard_dev_all_426_ad$sd_2
## z = 1.9272, n = 34, p-value = 0.05396
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## 131.0000000 4550.3333333 0.2335116
sd_sens_426_ad <- sens.slope(standard_dev_all_426_ad$sd_2)
print(sd_sens_426_ad)
##
## Sen's slope
##
## data: standard_dev_all_426_ad$sd_2
## z = 1.9272, n = 34, p-value = 0.05396
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.0008852415 0.0327236747
## sample estimates:
## Sen's slope
## 0.01613353
Oyler -> Morrisey 7/30/2003
snotel_457 <- SNOTEL_yampa_area %>%
filter(site_id == "457")
#str(snotel_457) # check the date, usually a character.
snotel_457$Date <- as.Date(snotel_457$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.
#THIS WILL CHANGE FOR EACH STATION
snotel_457_clean <- snotel_457 %>% # filter for the timeframe
filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
#filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers
addWaterYear() %>%
mutate(daymonth = format(as.Date(Date), "%d-%m")) %>%
na.omit()
#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))
snotel_457_clean <- snotel_457_clean %>%
group_by(waterYear)%>%
mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))
# Check for outliers
ggplot(snotel_457_clean, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
snotel_457_clean <- snotel_457_clean %>%
mutate(temp_diff = abs(temperature_min - temperature_max)) %>%
filter(temperature_mean > -40) %>%
filter(temperature_mean < 40)
ggplot(snotel_457_clean, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”
# filtering for temperature anomalies
#snotel_457_cull_count <- snotel_457_clean %>%
# filter(temperature_min > -40) %>%
# count(waterYear)
#snotel_457_cull_count
# filtering for too few observations in a year
snotel_457_cull_count_days <- snotel_457_clean %>%
group_by(waterYear) %>%
count(waterYear) %>%
filter(n < 350)
snotel_457_cull_count_days
## # A tibble: 6 x 2
## # Groups: waterYear [6]
## waterYear n
## <dbl> <int>
## 1 1985 1
## 2 1994 332
## 3 1996 349
## 4 1997 329
## 5 2003 346
## 6 2021 346
snotel_457_clean_culled <- snotel_457_clean %>%
filter(waterYear != "1985" & waterYear != "1986" & waterYear != "1994" & waterYear != "1996" & waterYear != "1997" & waterYear != "2003" & waterYear != "2021")# & waterYear != "1987" & waterYear != "1994" & waterYear != "2002")# & waterYear != "2002" & waterYear != "2016" & waterYear != "2022")# & waterYear != "2017") #%>%
#filter(temperature_mean > -49)
ggplot(snotel_457_clean_culled, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
Also 1986
ggplot(snotel_457_clean_culled, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
temp_457_xts <- xts(snotel_457_clean_culled$temperature_mean, order.by = snotel_457_clean_culled$Date)
dygraph(temp_457_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)")
#snotel_457_clean_culled <- snotel_457_clean_culled %>%
# filter(temperature_mean > -30)
#temp_457_xts <- xts(snotel_457_clean_culled$temperature_mean, order.by = snotel_457_clean_culled$Date)
#dygraph(temp_457_xts) %>%
# dyAxis("y", label = "Daily mean temperature (°C)")
snotel_457_adjusted <- snotel_457_clean_culled %>%
mutate(temp_ad = if_else(Date < "2003-07-30", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))
Non-corrected
#using the clean culled df:
#average water year temperature
yearly_wy_aver_457 <- snotel_457_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp = mean(temperature_mean))
#Average temperature by day for all water years:
daily_wy_aver_457 <- yearly_wy_aver_457 %>%
group_by(daymonth) %>%
mutate(aver_day_temp = mean(aver_ann_temp))
#average mean temperature by day for the period of record:
daily_wy_aver_457 <- daily_wy_aver_457 %>%
group_by(daymonth) %>%
mutate(all_ave_temp = mean(daily_wy_aver_457$aver_day_temp))
# try to show all years as means.
daily_wy_aver2_457 <-daily_wy_aver_457 %>%
group_by(waterDay) %>%
mutate(date_temp = mean(temperature_mean))
daily_wy_aver2_457$date_temp <- signif(daily_wy_aver2_457$date_temp,3) #reduce the sig figs
ggplot(daily_wy_aver2_457, aes(x = waterDay, y = date_temp))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
standard_dev_457 <- daily_wy_aver_457 %>%
group_by(waterYear) %>%
mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_457 <- standard_dev_457 %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_457 <- standard_dev_all_457 %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_457 %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')
| waterYear | sd_2 |
|---|---|
| 1987 | 8.892848 |
| 1988 | 9.953110 |
| 1989 | 9.757795 |
| 1990 | 9.176039 |
| 1991 | 9.421986 |
| 1992 | 8.507405 |
| 1993 | 8.878013 |
| 1995 | 8.358076 |
| 1998 | 9.150633 |
| 1999 | 8.440975 |
| 2000 | 8.736877 |
| 2001 | 9.645923 |
| 2002 | 10.040430 |
| 2004 | 8.051569 |
| 2005 | 7.863386 |
| 2006 | 8.768314 |
| 2007 | 9.036800 |
| 2008 | 8.900491 |
| 2009 | 7.976123 |
| 2010 | 8.677359 |
| 2011 | 8.610741 |
| 2012 | 8.805624 |
| 2013 | 9.212073 |
| 2014 | 8.350688 |
| 2015 | 7.879517 |
| 2016 | 8.589405 |
| 2017 | 8.326662 |
| 2018 | 8.323975 |
| 2019 | 8.654312 |
| 2020 | 9.071537 |
| 2022 | 8.555147 |
ggplot(standard_dev_all_457, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Non-corrected standard deviation of SNOTEL 457 average temperatures for water years 2005-2021
sd_mk_457 <- mk.test(standard_dev_all_457$sd_2)
print(sd_mk_457)
##
## Mann-Kendall trend test
##
## data: standard_dev_all_457$sd_2
## z = -2.3795, n = 31, p-value = 0.01734
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## -141.0000000 3461.6666667 -0.3032258
sd_sens_457 <- sens.slope(standard_dev_all_457$sd_2)
print(sd_sens_457)
##
## Sen's slope
##
## data: standard_dev_all_457$sd_2
## z = -2.3795, n = 31, p-value = 0.01734
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.049787462 -0.005081117
## sample estimates:
## Sen's slope
## -0.02638274
#using the clean culled df:
#average water year temperature
yearly_wy_aver_457_ad <- snotel_457_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_457_ad <- yearly_wy_aver_457_ad %>%
group_by(daymonth) %>%
mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_457_ad <- daily_wy_aver_457_ad %>%
group_by(daymonth) %>%
mutate(all_ave_temp_ad = mean(daily_wy_aver_457_ad$aver_day_temp_ad))
# try to show all years as means.
daily_wy_aver2_457_ad <-daily_wy_aver_457_ad %>%
group_by(waterDay) %>%
mutate(date_temp_ad = mean(temp_ad))
daily_wy_aver2_457_ad$date_temp_ad <- signif(daily_wy_aver2_457_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_457_ad, aes(x = waterDay, y = date_temp_ad))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
standard_dev_457_ad <- daily_wy_aver_457_ad %>%
group_by(waterYear) %>%
#filter(waterYear >= 1987 & waterYear <= 2021) %>%
mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_457_ad <- standard_dev_457_ad %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_457_ad <- standard_dev_all_457_ad %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_457_ad %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')
| waterYear | sd_2 |
|---|---|
| 1987 | 8.200529 |
| 1988 | 9.184255 |
| 1989 | 9.050691 |
| 1990 | 8.442843 |
| 1991 | 8.739959 |
| 1992 | 7.843690 |
| 1993 | 8.207632 |
| 1995 | 7.712921 |
| 1998 | 8.438042 |
| 1999 | 7.811811 |
| 2000 | 8.022352 |
| 2001 | 8.894507 |
| 2002 | 9.283490 |
| 2004 | 8.051478 |
| 2005 | 7.863308 |
| 2006 | 8.768371 |
| 2007 | 9.036790 |
| 2008 | 8.900641 |
| 2009 | 7.976349 |
| 2010 | 8.677249 |
| 2011 | 8.611094 |
| 2012 | 8.805285 |
| 2013 | 9.211821 |
| 2014 | 8.350803 |
| 2015 | 7.879602 |
| 2016 | 8.589409 |
| 2017 | 8.326313 |
| 2018 | 8.323899 |
| 2019 | 8.654521 |
| 2020 | 9.071250 |
| 2022 | 8.555395 |
ggplot(standard_dev_all_457_ad, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Morrisey-corrected standard deviation of SNOTEL 457 average temperatures for water years 1986-2021
sd_mk_457_ad <- mk.test(standard_dev_all_457_ad$sd_2)
print(sd_mk_457_ad)
##
## Mann-Kendall trend test
##
## data: standard_dev_all_457_ad$sd_2
## z = 0.44191, n = 31, p-value = 0.6586
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## 2.700000e+01 3.461667e+03 5.806452e-02
sd_sens_457_ad <- sens.slope(standard_dev_all_457_ad$sd_2)
print(sd_sens_457_ad)
##
## Sen's slope
##
## data: standard_dev_all_457_ad$sd_2
## z = 0.44191, n = 31, p-value = 0.6586
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.01772587 0.02844930
## sample estimates:
## Sen's slope
## 0.004837867
Oyler -> Morrisey 8/7/2006
snotel_467 <- SNOTEL_yampa_area %>%
filter(site_id == "467")
#str(snotel_467) # check the date, usually a character.
snotel_467$Date <- as.Date(snotel_467$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.
#THIS WILL CHANGE FOR EACH STATION
snotel_467_clean <- snotel_467 %>% # filter for the timeframe
filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
#filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers
addWaterYear() %>%
mutate(daymonth = format(as.Date(Date), "%d-%m")) %>%
na.omit()
#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))
snotel_467_clean <- snotel_467_clean %>%
group_by(waterYear)%>%
mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))
# Check for outliers
ggplot(snotel_467_clean, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
snotel_467_clean <- snotel_467_clean %>%
mutate(temp_diff = abs(temperature_min - temperature_max)) %>%
filter(temperature_mean > -40) %>%
filter(temperature_mean < 40)
ggplot(snotel_467_clean, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”
# filtering for temperature anomalies
#snotel_467_cull_count <- snotel_467_clean %>%
# filter(temperature_min > -40) %>%
# count(waterYear)
#snotel_467_cull_count
# filtering for too few observations in a year
snotel_467_cull_count_days <- snotel_467_clean %>%
group_by(waterYear) %>%
count(waterYear) %>%
filter(n < 350)
snotel_467_cull_count_days
## # A tibble: 6 x 2
## # Groups: waterYear [6]
## waterYear n
## <dbl> <int>
## 1 1985 1
## 2 1994 326
## 3 1998 342
## 4 1999 281
## 5 2001 339
## 6 2021 344
snotel_467_clean_culled <- snotel_467_clean %>%
filter(waterYear != "1985" & waterYear != "1986" & waterYear != "1994" & waterYear != "1998" & waterYear != "1999" & waterYear != "2001" & waterYear != "2021")# & waterYear != "2021")# & waterYear != "1987" & waterYear != "1994" & waterYear != "2002")# & waterYear != "2002" & waterYear != "2016" & waterYear != "2022")# & waterYear != "2017") #%>%
#filter(temperature_mean > -49)
ggplot(snotel_467_clean_culled, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
Also 1986
ggplot(snotel_467_clean_culled, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
temp_467_xts <- xts(snotel_467_clean_culled$temperature_mean, order.by = snotel_467_clean_culled$Date)
dygraph(temp_467_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)")
#snotel_467_clean_culled <- snotel_467_clean_culled %>%
# filter(temperature_mean > -30)
#temp_467_xts <- xts(snotel_467_clean_culled$temperature_mean, order.by = snotel_467_clean_culled$Date)
#dygraph(temp_467_xts) %>%
# dyAxis("y", label = "Daily mean temperature (°C)")
snotel_467_adjusted <- snotel_467_clean_culled %>%
mutate(temp_ad = if_else(Date < "2006-08-07", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))
Non-corrected
#using the clean culled df:
#average water year temperature
yearly_wy_aver_467 <- snotel_467_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp = mean(temperature_mean))
#Average temperature by day for all water years:
daily_wy_aver_467 <- yearly_wy_aver_467 %>%
group_by(daymonth) %>%
mutate(aver_day_temp = mean(aver_ann_temp))
#average mean temperature by day for the period of record:
daily_wy_aver_467 <- daily_wy_aver_467 %>%
group_by(daymonth) %>%
mutate(all_ave_temp = mean(daily_wy_aver_467$aver_day_temp))
# try to show all years as means.
daily_wy_aver2_467 <-daily_wy_aver_467 %>%
group_by(waterDay) %>%
mutate(date_temp = mean(temperature_mean))
daily_wy_aver2_467$date_temp <- signif(daily_wy_aver2_467$date_temp,3) #reduce the sig figs
ggplot(daily_wy_aver2_467, aes(x = waterDay, y = date_temp))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
standard_dev_467 <- daily_wy_aver_467 %>%
group_by(waterYear) %>%
mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_467 <- standard_dev_467 %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_467 <- standard_dev_all_467 %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_467 %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')
| waterYear | sd_2 |
|---|---|
| 1987 | 9.103406 |
| 1988 | 10.510726 |
| 1989 | 10.049523 |
| 1990 | 9.547427 |
| 1991 | 9.692602 |
| 1992 | 8.717383 |
| 1993 | 9.114798 |
| 1995 | 8.625842 |
| 1996 | 8.985964 |
| 1997 | 9.119643 |
| 2000 | 8.850761 |
| 2002 | 10.120172 |
| 2003 | 9.161089 |
| 2004 | 8.589463 |
| 2005 | 8.525521 |
| 2006 | 9.564659 |
| 2007 | 9.088368 |
| 2008 | 9.090316 |
| 2009 | 8.102303 |
| 2010 | 8.886760 |
| 2011 | 8.734057 |
| 2012 | 8.902720 |
| 2013 | 9.375328 |
| 2014 | 8.477094 |
| 2015 | 8.044768 |
| 2016 | 8.660871 |
| 2017 | 8.514442 |
| 2018 | 8.465435 |
| 2019 | 8.883368 |
| 2020 | 9.345838 |
| 2022 | 8.902383 |
ggplot(standard_dev_all_467, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Non-corrected standard deviation of SNOTEL 467 average temperatures for water years 2005-2021
sd_mk_467 <- mk.test(standard_dev_all_467$sd_2)
print(sd_mk_467)
##
## Mann-Kendall trend test
##
## data: standard_dev_all_467$sd_2
## z = -2.7194, n = 31, p-value = 0.00654
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## -161.0000000 3461.6666667 -0.3462366
sd_sens_467 <- sens.slope(standard_dev_all_467$sd_2)
print(sd_sens_467)
##
## Sen's slope
##
## data: standard_dev_all_467$sd_2
## z = -2.7194, n = 31, p-value = 0.00654
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.053188331 -0.007753413
## sample estimates:
## Sen's slope
## -0.02739628
#using the clean culled df:
#average water year temperature
yearly_wy_aver_467_ad <- snotel_467_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_467_ad <- yearly_wy_aver_467_ad %>%
group_by(daymonth) %>%
mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_467_ad <- daily_wy_aver_467_ad %>%
group_by(daymonth) %>%
mutate(all_ave_temp_ad = mean(daily_wy_aver_467_ad$aver_day_temp_ad))
# try to show all years as means.
daily_wy_aver2_467_ad <-daily_wy_aver_467_ad %>%
group_by(waterDay) %>%
mutate(date_temp_ad = mean(temp_ad))
daily_wy_aver2_467_ad$date_temp_ad <- signif(daily_wy_aver2_467_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_467_ad, aes(x = waterDay, y = date_temp_ad))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
standard_dev_467_ad <- daily_wy_aver_467_ad %>%
group_by(waterYear) %>%
#filter(waterYear >= 1987 & waterYear <= 2021) %>%
mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_467_ad <- standard_dev_467_ad %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_467_ad <- standard_dev_all_467_ad %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_467_ad %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')
| waterYear | sd_2 |
|---|---|
| 1987 | 8.416565 |
| 1988 | 9.721370 |
| 1989 | 9.348058 |
| 1990 | 8.816292 |
| 1991 | 9.016762 |
| 1992 | 8.057564 |
| 1993 | 8.438420 |
| 1995 | 7.963510 |
| 1996 | 8.304891 |
| 1997 | 8.454695 |
| 2000 | 8.126220 |
| 2002 | 9.357756 |
| 2003 | 8.423881 |
| 2004 | 7.952025 |
| 2005 | 7.847462 |
| 2006 | 8.843490 |
| 2007 | 9.088080 |
| 2008 | 9.090597 |
| 2009 | 8.102286 |
| 2010 | 8.886876 |
| 2011 | 8.734381 |
| 2012 | 8.902800 |
| 2013 | 9.375482 |
| 2014 | 8.477382 |
| 2015 | 8.045353 |
| 2016 | 8.661195 |
| 2017 | 8.514535 |
| 2018 | 8.465471 |
| 2019 | 8.883707 |
| 2020 | 9.345830 |
| 2022 | 8.902408 |
ggplot(standard_dev_all_467_ad, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Morrisey-corrected standard deviation of SNOTEL 467 average temperatures for water years 1986-2021
sd_mk_467_ad <- mk.test(standard_dev_all_467_ad$sd_2)
print(sd_mk_467_ad)
##
## Mann-Kendall trend test
##
## data: standard_dev_all_467_ad$sd_2
## z = 0.57788, n = 31, p-value = 0.5633
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## 3.500000e+01 3.461667e+03 7.526882e-02
sd_sens_467_ad <- sens.slope(standard_dev_all_467_ad$sd_2)
print(sd_sens_467_ad)
##
## Sen's slope
##
## data: standard_dev_all_467_ad$sd_2
## z = 0.57788, n = 31, p-value = 0.5633
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.01785969 0.02658485
## sample estimates:
## Sen's slope
## 0.003440178
Morrisey 5/22/2006
snotel_607 <- SNOTEL_yampa_area %>%
filter(site_id == "607")
#str(snotel_607) # check the date, usually a character.
snotel_607$Date <- as.Date(snotel_607$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.
#THIS WILL CHANGE FOR EACH STATION
snotel_607_clean <- snotel_607 %>% # filter for the timeframe
filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
#filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers
addWaterYear() %>%
mutate(daymonth = format(as.Date(Date), "%d-%m")) %>%
na.omit()
#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))
snotel_607_clean <- snotel_607_clean %>%
group_by(waterYear)%>%
mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))
# Check for outliers
ggplot(snotel_607_clean, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
snotel_607_clean <- snotel_607_clean %>%
mutate(temp_diff = abs(temperature_min - temperature_max)) %>%
filter(temperature_mean > -40) %>%
filter(temperature_mean < 40)
ggplot(snotel_607_clean, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”
# filtering for temperature anomalies
#snotel_607_cull_count <- snotel_607_clean %>%
# filter(temperature_min > -40) %>%
# count(waterYear)
#snotel_607_cull_count
# filtering for too few observations in a year
snotel_607_cull_count_days <- snotel_607_clean %>%
group_by(waterYear) %>%
count(waterYear) %>%
filter(n < 350)
snotel_607_cull_count_days
## # A tibble: 4 x 2
## # Groups: waterYear [4]
## waterYear n
## <dbl> <int>
## 1 1985 1
## 2 1994 333
## 3 2019 284
## 4 2020 328
snotel_607_clean_culled <- snotel_607_clean %>%
filter(waterYear != "1985" & waterYear != "1986" & waterYear != "1994" & waterYear != "2019" & waterYear != "2020")# & waterYear != "2001" & waterYear != "2021")# & waterYear != "2021")# & waterYear != "1987" & waterYear != "1994" & waterYear != "2002")# & waterYear != "2002" & waterYear != "2016" & waterYear != "2022")# & waterYear != "2017") #%>%
#filter(temperature_mean > -49)
ggplot(snotel_607_clean_culled, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
Also 1986
ggplot(snotel_607_clean_culled, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
temp_607_xts <- xts(snotel_607_clean_culled$temperature_mean, order.by = snotel_607_clean_culled$Date)
dygraph(temp_607_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)")
#snotel_607_clean_culled <- snotel_607_clean_culled %>%
# filter(temperature_mean > -30)
#temp_607_xts <- xts(snotel_607_clean_culled$temperature_mean, order.by = snotel_607_clean_culled$Date)
#dygraph(temp_607_xts) %>%
# dyAxis("y", label = "Daily mean temperature (°C)")
snotel_607_adjusted <- snotel_607_clean_culled %>%
mutate(temp_ad = if_else(Date < "2006-05-22", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))
Non-corrected
#using the clean culled df:
#average water year temperature
yearly_wy_aver_607 <- snotel_607_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp = mean(temperature_mean))
#Average temperature by day for all water years:
daily_wy_aver_607 <- yearly_wy_aver_607 %>%
group_by(daymonth) %>%
mutate(aver_day_temp = mean(aver_ann_temp))
#average mean temperature by day for the period of record:
daily_wy_aver_607 <- daily_wy_aver_607 %>%
group_by(daymonth) %>%
mutate(all_ave_temp = mean(daily_wy_aver_607$aver_day_temp))
# try to show all years as means.
daily_wy_aver2_607 <-daily_wy_aver_607 %>%
group_by(waterDay) %>%
mutate(date_temp = mean(temperature_mean))
daily_wy_aver2_607$date_temp <- signif(daily_wy_aver2_607$date_temp,3) #reduce the sig figs
ggplot(daily_wy_aver2_607, aes(x = waterDay, y = date_temp))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
standard_dev_607 <- daily_wy_aver_607 %>%
group_by(waterYear) %>%
mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_607 <- standard_dev_607 %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_607 <- standard_dev_all_607 %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_607 %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')
| waterYear | sd_2 |
|---|---|
| 1987 | 8.658921 |
| 1988 | 9.691007 |
| 1989 | 9.677093 |
| 1990 | 9.127534 |
| 1991 | 9.267600 |
| 1992 | 8.495617 |
| 1993 | 8.542784 |
| 1995 | 8.427463 |
| 1996 | 8.936524 |
| 1997 | 8.897439 |
| 1998 | 9.097517 |
| 1999 | 8.407095 |
| 2000 | 8.685669 |
| 2001 | 9.559462 |
| 2002 | 10.100422 |
| 2003 | 8.935295 |
| 2004 | 8.581025 |
| 2005 | 8.571462 |
| 2006 | 9.563477 |
| 2007 | 9.027055 |
| 2008 | 8.751128 |
| 2009 | 7.935985 |
| 2010 | 9.022507 |
| 2011 | 8.612122 |
| 2012 | 8.812590 |
| 2013 | 9.410619 |
| 2014 | 8.353556 |
| 2015 | 7.907493 |
| 2016 | 8.522596 |
| 2017 | 8.366193 |
| 2018 | 8.240463 |
| 2021 | 8.972977 |
| 2022 | 8.767136 |
ggplot(standard_dev_all_607, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Non-corrected standard deviation of SNOTEL 607 average temperatures for water years 2005-2021
sd_mk_607 <- mk.test(standard_dev_all_607$sd_2)
print(sd_mk_607)
##
## Mann-Kendall trend test
##
## data: standard_dev_all_607$sd_2
## z = -2.0298, n = 33, p-value = 0.04238
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## -132.000 4165.333 -0.250
sd_sens_607 <- sens.slope(standard_dev_all_607$sd_2)
print(sd_sens_607)
##
## Sen's slope
##
## data: standard_dev_all_607$sd_2
## z = -2.0298, n = 33, p-value = 0.04238
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.036056261 -0.001515865
## sample estimates:
## Sen's slope
## -0.01788946
#using the clean culled df:
#average water year temperature
yearly_wy_aver_607_ad <- snotel_607_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_607_ad <- yearly_wy_aver_607_ad %>%
group_by(daymonth) %>%
mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_607_ad <- daily_wy_aver_607_ad %>%
group_by(daymonth) %>%
mutate(all_ave_temp_ad = mean(daily_wy_aver_607_ad$aver_day_temp_ad))
# try to show all years as means.
daily_wy_aver2_607_ad <-daily_wy_aver_607_ad %>%
group_by(waterDay) %>%
mutate(date_temp_ad = mean(temp_ad))
daily_wy_aver2_607_ad$date_temp_ad <- signif(daily_wy_aver2_607_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_607_ad, aes(x = waterDay, y = date_temp_ad))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
standard_dev_607_ad <- daily_wy_aver_607_ad %>%
group_by(waterYear) %>%
#filter(waterYear >= 1987 & waterYear <= 2021) %>%
mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_607_ad <- standard_dev_607_ad %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_607_ad <- standard_dev_all_607_ad %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_607_ad %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')
| waterYear | sd_2 |
|---|---|
| 1987 | 8.065805 |
| 1988 | 9.041423 |
| 1989 | 9.118850 |
| 1990 | 8.514101 |
| 1991 | 8.699806 |
| 1992 | 7.927056 |
| 1993 | 7.982533 |
| 1995 | 7.846909 |
| 1996 | 8.320426 |
| 1997 | 8.291259 |
| 1998 | 8.445612 |
| 1999 | 7.829298 |
| 2000 | 8.031514 |
| 2001 | 8.877490 |
| 2002 | 9.409575 |
| 2003 | 8.273822 |
| 2004 | 8.002022 |
| 2005 | 7.966244 |
| 2006 | 8.852307 |
| 2007 | 9.027252 |
| 2008 | 8.750916 |
| 2009 | 7.936108 |
| 2010 | 9.022503 |
| 2011 | 8.612197 |
| 2012 | 8.812993 |
| 2013 | 9.410822 |
| 2014 | 8.353805 |
| 2015 | 7.907862 |
| 2016 | 8.523025 |
| 2017 | 8.366000 |
| 2018 | 8.241006 |
| 2021 | 8.973005 |
| 2022 | 8.767834 |
ggplot(standard_dev_all_607_ad, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Morrisey-corrected standard deviation of SNOTEL 607 average temperatures for water years 1986-2021
sd_mk_607_ad <- mk.test(standard_dev_all_607_ad$sd_2)
print(sd_mk_607_ad)
##
## Mann-Kendall trend test
##
## data: standard_dev_all_607_ad$sd_2
## z = 0.63527, n = 33, p-value = 0.5253
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## 4.200000e+01 4.165333e+03 7.954545e-02
sd_sens_607_ad <- sens.slope(standard_dev_all_607_ad$sd_2)
print(sd_sens_607_ad)
##
## Sen's slope
##
## data: standard_dev_all_607_ad$sd_2
## z = 0.63527, n = 33, p-value = 0.5253
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.008825448 0.026678756
## sample estimates:
## Sen's slope
## 0.006477733
Morrisey 8/7/2006
snotel_709 <- SNOTEL_yampa_area %>%
filter(site_id == "709")
#str(snotel_709) # check the date, usually a character.
snotel_709$Date <- as.Date(snotel_709$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.
#THIS WILL CHANGE FOR EACH STATION
snotel_709_clean <- snotel_709 %>% # filter for the timeframe
filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
#filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers
addWaterYear() %>%
mutate(daymonth = format(as.Date(Date), "%d-%m")) %>%
na.omit()
#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))
snotel_709_clean <- snotel_709_clean %>%
group_by(waterYear)%>%
mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))
# Check for outliers
ggplot(snotel_709_clean, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
snotel_709_clean <- snotel_709_clean %>%
mutate(temp_diff = abs(temperature_min - temperature_max)) %>%
filter(temperature_mean > -40) %>%
filter(temperature_mean < 40)
ggplot(snotel_709_clean, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”
# filtering for temperature anomalies
#snotel_709_cull_count <- snotel_709_clean %>%
# filter(temperature_min > -40) %>%
# count(waterYear)
#snotel_709_cull_count
# filtering for too few observations in a year
snotel_709_cull_count_days <- snotel_709_clean %>%
group_by(waterYear) %>%
count(waterYear) %>%
filter(n < 350)
snotel_709_cull_count_days
## # A tibble: 14 x 2
## # Groups: waterYear [14]
## waterYear n
## <dbl> <int>
## 1 1992 318
## 2 1993 293
## 3 1994 332
## 4 1995 169
## 5 2006 348
## 6 2007 311
## 7 2011 292
## 8 2012 308
## 9 2013 300
## 10 2014 283
## 11 2015 44
## 12 2016 18
## 13 2017 34
## 14 2018 147
snotel_709_clean_culled <- snotel_709_clean %>%
filter(waterYear > "1995" & waterYear != "2006" & waterYear != "2007" & waterYear != "2011" & waterYear != "2012" & waterYear != "2013" & waterYear != "2014" & waterYear != "2015" & waterYear != "2016" & waterYear != "2017" & waterYear != "2018")# & waterYear != "2002")# & waterYear != "2002" & waterYear != "2016" & waterYear != "2022")# & waterYear != "2017") #%>%
#filter(temperature_mean > -49)
ggplot(snotel_709_clean_culled, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
ggplot(snotel_709_clean_culled, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
temp_709_xts <- xts(snotel_709_clean_culled$temperature_mean, order.by = snotel_709_clean_culled$Date)
dygraph(temp_709_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)")
#snotel_709_clean_culled <- snotel_709_clean_culled %>%
# filter(temperature_mean > -30)
#temp_709_xts <- xts(snotel_709_clean_culled$temperature_mean, order.by = snotel_709_clean_culled$Date)
#dygraph(temp_709_xts) %>%
# dyAxis("y", label = "Daily mean temperature (°C)")
snotel_709_adjusted <- snotel_709_clean_culled %>%
mutate(temp_ad = if_else(Date < "2006-08-07", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))
Non-corrected
#using the clean culled df:
#average water year temperature
yearly_wy_aver_709 <- snotel_709_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp = mean(temperature_mean))
#Average temperature by day for all water years:
daily_wy_aver_709 <- yearly_wy_aver_709 %>%
group_by(daymonth) %>%
mutate(aver_day_temp = mean(aver_ann_temp))
#average mean temperature by day for the period of record:
daily_wy_aver_709 <- daily_wy_aver_709 %>%
group_by(daymonth) %>%
mutate(all_ave_temp = mean(daily_wy_aver_709$aver_day_temp))
# try to show all years as means.
daily_wy_aver2_709 <-daily_wy_aver_709 %>%
group_by(waterDay) %>%
mutate(date_temp = mean(temperature_mean))
daily_wy_aver2_709$date_temp <- signif(daily_wy_aver2_709$date_temp,3) #reduce the sig figs
ggplot(daily_wy_aver2_709, aes(x = waterDay, y = date_temp))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
standard_dev_709 <- daily_wy_aver_709 %>%
group_by(waterYear) %>%
mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_709 <- standard_dev_709 %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_709 <- standard_dev_all_709 %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_709 %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')
| waterYear | sd_2 |
|---|---|
| 1996 | 8.845436 |
| 1997 | 8.864608 |
| 1998 | 9.079136 |
| 1999 | 8.309577 |
| 2000 | 8.748503 |
| 2001 | 9.426260 |
| 2002 | 9.843776 |
| 2003 | 9.011443 |
| 2004 | 8.578795 |
| 2005 | 8.457080 |
| 2008 | 9.038776 |
| 2009 | 8.123905 |
| 2010 | 8.958458 |
| 2019 | 8.740431 |
| 2020 | 9.196584 |
| 2021 | 9.081772 |
| 2022 | 8.759943 |
ggplot(standard_dev_all_709, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Non-corrected standard deviation of SNOTEL 709 average temperatures for water years 2005-2021
sd_mk_709 <- mk.test(standard_dev_all_709$sd_2)
print(sd_mk_709)
##
## Mann-Kendall trend test
##
## data: standard_dev_all_709$sd_2
## z = 0.041193, n = 17, p-value = 0.9671
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## 2.00000000 589.33333333 0.01470588
sd_sens_709 <- sens.slope(standard_dev_all_709$sd_2)
print(sd_sens_709)
##
## Sen's slope
##
## data: standard_dev_all_709$sd_2
## z = 0.041193, n = 17, p-value = 0.9671
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.04647212 0.03945218
## sample estimates:
## Sen's slope
## 0.0005780099
#using the clean culled df:
#average water year temperature
yearly_wy_aver_709_ad <- snotel_709_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_709_ad <- yearly_wy_aver_709_ad %>%
group_by(daymonth) %>%
mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_709_ad <- daily_wy_aver_709_ad %>%
group_by(daymonth) %>%
mutate(all_ave_temp_ad = mean(daily_wy_aver_709_ad$aver_day_temp_ad))
# try to show all years as means.
daily_wy_aver2_709_ad <-daily_wy_aver_709_ad %>%
group_by(waterDay) %>%
mutate(date_temp_ad = mean(temp_ad))
daily_wy_aver2_709_ad$date_temp_ad <- signif(daily_wy_aver2_709_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_709_ad, aes(x = waterDay, y = date_temp_ad))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
standard_dev_709_ad <- daily_wy_aver_709_ad %>%
group_by(waterYear) %>%
#filter(waterYear >= 1987 & waterYear <= 2021) %>%
mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_709_ad <- standard_dev_709_ad %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_709_ad <- standard_dev_all_709_ad %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_709_ad %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')
| waterYear | sd_2 |
|---|---|
| 1996 | 8.213076 |
| 1997 | 8.255603 |
| 1998 | 8.412792 |
| 1999 | 7.726856 |
| 2000 | 8.073670 |
| 2001 | 8.727586 |
| 2002 | 9.142221 |
| 2003 | 8.322566 |
| 2004 | 7.979937 |
| 2005 | 7.819602 |
| 2008 | 9.038980 |
| 2009 | 8.124255 |
| 2010 | 8.958543 |
| 2019 | 8.740727 |
| 2020 | 9.196244 |
| 2021 | 9.081146 |
| 2022 | 8.760212 |
ggplot(standard_dev_all_709_ad, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Morrisey-corrected standard deviation of SNOTEL 709 average temperatures for water years 1986-2021
sd_mk_709_ad <- mk.test(standard_dev_all_709_ad$sd_2)
print(sd_mk_709_ad)
##
## Mann-Kendall trend test
##
## data: standard_dev_all_709_ad$sd_2
## z = 1.9361, n = 17, p-value = 0.05286
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## 48.0000000 589.3333333 0.3529412
sd_sens_709_ad <- sens.slope(standard_dev_all_709_ad$sd_2)
print(sd_sens_709_ad)
##
## Sen's slope
##
## data: standard_dev_all_709_ad$sd_2
## z = 1.9361, n = 17, p-value = 0.05286
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.006786111 0.102902003
## sample estimates:
## Sen's slope
## 0.05101403
Oyler -> Morrisey 8/7/2006
snotel_717 <- SNOTEL_yampa_area %>%
filter(site_id == "717")
#str(snotel_717) # check the date, usually a character.
snotel_717$Date <- as.Date(snotel_717$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.
#THIS WILL CHANGE FOR EACH STATION
snotel_717_clean <- snotel_717 %>% # filter for the timeframe
filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
#filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers
addWaterYear() %>%
mutate(daymonth = format(as.Date(Date), "%d-%m")) %>%
na.omit()
#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))
snotel_717_clean <- snotel_717_clean %>%
group_by(waterYear)%>%
mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))
# Check for outliers
ggplot(snotel_717_clean, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
snotel_717_clean <- snotel_717_clean %>%
mutate(temp_diff = abs(temperature_min - temperature_max)) %>%
filter(temperature_mean > -40) %>%
filter(temperature_mean < 40)
ggplot(snotel_717_clean, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”
# filtering for temperature anomalies
#snotel_717_cull_count <- snotel_717_clean %>%
# filter(temperature_min > -40) %>%
# count(waterYear)
#snotel_717_cull_count
# filtering for too few observations in a year
snotel_717_cull_count_days <- snotel_717_clean %>%
group_by(waterYear) %>%
count(waterYear) %>%
filter(n < 350)
snotel_717_cull_count_days
## # A tibble: 2 x 2
## # Groups: waterYear [2]
## waterYear n
## <dbl> <int>
## 1 1992 311
## 2 1993 330
snotel_717_clean_culled <- snotel_717_clean %>%
filter(waterYear != "1992" & waterYear != "1993")# & waterYear != "2007" & waterYear != "2011" & waterYear != "2012" & waterYear != "2013" & waterYear != "2014" & waterYear != "2015" & waterYear != "2016" & waterYear != "2017" & waterYear != "2018")# & waterYear != "2002")# & waterYear != "2002" & waterYear != "2016" & waterYear != "2022")# & waterYear != "2017") #%>%
#filter(temperature_mean > -49)
ggplot(snotel_717_clean_culled, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
ggplot(snotel_717_clean_culled, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
temp_717_xts <- xts(snotel_717_clean_culled$temperature_mean, order.by = snotel_717_clean_culled$Date)
dygraph(temp_717_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)")
#snotel_717_clean_culled <- snotel_717_clean_culled %>%
# filter(temperature_mean > -30)
#temp_717_xts <- xts(snotel_717_clean_culled$temperature_mean, order.by = snotel_717_clean_culled$Date)
#dygraph(temp_717_xts) %>%
# dyAxis("y", label = "Daily mean temperature (°C)")
snotel_717_adjusted <- snotel_717_clean_culled %>%
mutate(temp_ad = if_else(Date < "2006-08-07", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))
Non-corrected
#using the clean culled df:
#average water year temperature
yearly_wy_aver_717 <- snotel_717_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp = mean(temperature_mean))
#Average temperature by day for all water years:
daily_wy_aver_717 <- yearly_wy_aver_717 %>%
group_by(daymonth) %>%
mutate(aver_day_temp = mean(aver_ann_temp))
#average mean temperature by day for the period of record:
daily_wy_aver_717 <- daily_wy_aver_717 %>%
group_by(daymonth) %>%
mutate(all_ave_temp = mean(daily_wy_aver_717$aver_day_temp))
# try to show all years as means.
daily_wy_aver2_717 <-daily_wy_aver_717 %>%
group_by(waterDay) %>%
mutate(date_temp = mean(temperature_mean))
daily_wy_aver2_717$date_temp <- signif(daily_wy_aver2_717$date_temp,3) #reduce the sig figs
ggplot(daily_wy_aver2_717, aes(x = waterDay, y = date_temp))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
standard_dev_717 <- daily_wy_aver_717 %>%
group_by(waterYear) %>%
mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_717 <- standard_dev_717 %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_717 <- standard_dev_all_717 %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_717 %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')
| waterYear | sd_2 |
|---|---|
| 1987 | 8.803468 |
| 1988 | 9.707652 |
| 1989 | 9.270222 |
| 1990 | 9.163954 |
| 1991 | 9.043756 |
| 1994 | 9.475534 |
| 1995 | 8.216129 |
| 1996 | 8.766674 |
| 1997 | 8.761084 |
| 1998 | 8.941067 |
| 1999 | 8.106802 |
| 2000 | 8.669178 |
| 2001 | 9.209650 |
| 2002 | 9.676307 |
| 2003 | 9.008971 |
| 2004 | 8.362946 |
| 2005 | 8.373378 |
| 2006 | 9.248617 |
| 2007 | 8.800246 |
| 2008 | 8.921082 |
| 2009 | 7.998549 |
| 2010 | 8.681956 |
| 2011 | 8.605749 |
| 2012 | 8.582569 |
| 2013 | 9.046302 |
| 2014 | 8.329385 |
| 2015 | 7.737981 |
| 2016 | 8.402920 |
| 2017 | 8.250374 |
| 2018 | 8.202730 |
| 2019 | 8.544093 |
| 2020 | 8.992894 |
| 2021 | 8.957109 |
| 2022 | 8.519914 |
ggplot(standard_dev_all_717, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Non-corrected standard deviation of SNOTEL 717 average temperatures for water years 2005-2021
sd_mk_717 <- mk.test(standard_dev_all_717$sd_2)
print(sd_mk_717)
##
## Mann-Kendall trend test
##
## data: standard_dev_all_717$sd_2
## z = -2.6684, n = 34, p-value = 0.007621
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## -181.0000000 4550.3333333 -0.3226381
sd_sens_717 <- sens.slope(standard_dev_all_717$sd_2)
print(sd_sens_717)
##
## Sen's slope
##
## data: standard_dev_all_717$sd_2
## z = -2.6684, n = 34, p-value = 0.007621
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.039224947 -0.006265457
## sample estimates:
## Sen's slope
## -0.02259117
#using the clean culled df:
#average water year temperature
yearly_wy_aver_717_ad <- snotel_717_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_717_ad <- yearly_wy_aver_717_ad %>%
group_by(daymonth) %>%
mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_717_ad <- daily_wy_aver_717_ad %>%
group_by(daymonth) %>%
mutate(all_ave_temp_ad = mean(daily_wy_aver_717_ad$aver_day_temp_ad))
# try to show all years as means.
daily_wy_aver2_717_ad <-daily_wy_aver_717_ad %>%
group_by(waterDay) %>%
mutate(date_temp_ad = mean(temp_ad))
daily_wy_aver2_717_ad$date_temp_ad <- signif(daily_wy_aver2_717_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_717_ad, aes(x = waterDay, y = date_temp_ad))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
standard_dev_717_ad <- daily_wy_aver_717_ad %>%
group_by(waterYear) %>%
#filter(waterYear >= 1987 & waterYear <= 2021) %>%
mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_717_ad <- standard_dev_717_ad %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_717_ad <- standard_dev_all_717_ad %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_717_ad %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')
| waterYear | sd_2 |
|---|---|
| 1987 | 8.203857 |
| 1988 | 9.063074 |
| 1989 | 8.683957 |
| 1990 | 8.529404 |
| 1991 | 8.477488 |
| 1994 | 8.832494 |
| 1995 | 7.643847 |
| 1996 | 8.162790 |
| 1997 | 8.181297 |
| 1998 | 8.315921 |
| 1999 | 7.566445 |
| 2000 | 8.038773 |
| 2001 | 8.563752 |
| 2002 | 9.013927 |
| 2003 | 8.355528 |
| 2004 | 7.809021 |
| 2005 | 7.779504 |
| 2006 | 8.600879 |
| 2007 | 8.799659 |
| 2008 | 8.920535 |
| 2009 | 7.997923 |
| 2010 | 8.681418 |
| 2011 | 8.605311 |
| 2012 | 8.581992 |
| 2013 | 9.045623 |
| 2014 | 8.328989 |
| 2015 | 7.737670 |
| 2016 | 8.402569 |
| 2017 | 8.249042 |
| 2018 | 8.202635 |
| 2019 | 8.543763 |
| 2020 | 8.992050 |
| 2021 | 8.956530 |
| 2022 | 8.520306 |
ggplot(standard_dev_all_717_ad, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Morrisey-corrected standard deviation of SNOTEL 717 average temperatures for water years 1986-2021
sd_mk_717_ad <- mk.test(standard_dev_all_717_ad$sd_2)
print(sd_mk_717_ad)
##
## Mann-Kendall trend test
##
## data: standard_dev_all_717_ad$sd_2
## z = 0.56333, n = 34, p-value = 0.5732
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## 3.900000e+01 4.550333e+03 6.951872e-02
sd_sens_717_ad <- sens.slope(standard_dev_all_717_ad$sd_2)
print(sd_sens_717_ad)
##
## Sen's slope
##
## data: standard_dev_all_717_ad$sd_2
## z = 0.56333, n = 34, p-value = 0.5732
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.01114957 0.02335425
## sample estimates:
## Sen's slope
## 0.004691148
Original 8/18/2004
snotel_825 <- SNOTEL_yampa_area %>%
filter(site_id == "825")
#str(snotel_825) # check the date, usually a character.
snotel_825$Date <- as.Date(snotel_825$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.
#THIS WILL CHANGE FOR EACH STATION
snotel_825_clean <- snotel_825 %>% # filter for the timeframe
filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
#filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers
addWaterYear() %>%
mutate(daymonth = format(as.Date(Date), "%d-%m")) %>%
na.omit()
#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))
snotel_825_clean <- snotel_825_clean %>%
group_by(waterYear)%>%
mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))
# Check for outliers
ggplot(snotel_825_clean, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
snotel_825_clean <- snotel_825_clean %>%
mutate(temp_diff = abs(temperature_min - temperature_max)) %>%
filter(temperature_mean > -40) %>%
filter(temperature_mean < 30)
ggplot(snotel_825_clean, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”
# filtering for temperature anomalies
#snotel_825_cull_count <- snotel_825_clean %>%
# filter(temperature_min > -40) %>%
# count(waterYear)
#snotel_825_cull_count
# filtering for too few observations in a year
snotel_825_cull_count_days <- snotel_825_clean %>%
group_by(waterYear) %>%
count(waterYear) %>%
filter(n < 350)
snotel_825_cull_count_days
## # A tibble: 9 x 2
## # Groups: waterYear [9]
## waterYear n
## <dbl> <int>
## 1 1985 1
## 2 1987 348
## 3 1994 301
## 4 1995 347
## 5 1997 322
## 6 1998 343
## 7 1999 348
## 8 2002 327
## 9 2008 325
snotel_825_clean_culled <- snotel_825_clean %>%
filter(waterYear > "1987" & waterYear != "1994" & waterYear != "1995" & waterYear != "1997" & waterYear != "1998" & waterYear != "1999" & waterYear != "2002" & waterYear != "2008")# & waterYear != "2017" & waterYear != "2018")# & waterYear != "2002")# & waterYear != "2002" & waterYear != "2016" & waterYear != "2022")# & waterYear != "2017") #%>%
#filter(temperature_mean > -49)
ggplot(snotel_825_clean_culled, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
ggplot(snotel_825_clean_culled, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
temp_825_xts <- xts(snotel_825_clean_culled$temperature_mean, order.by = snotel_825_clean_culled$Date)
dygraph(temp_825_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)")
#snotel_825_clean_culled <- snotel_825_clean_culled %>%
# filter(temperature_mean > -30)
#temp_825_xts <- xts(snotel_825_clean_culled$temperature_mean, order.by = snotel_825_clean_culled$Date)
#dygraph(temp_825_xts) %>%
# dyAxis("y", label = "Daily mean temperature (°C)")
snotel_825_adjusted <- snotel_825_clean_culled %>%
mutate(temp_ad = if_else(Date < "2004-08-18", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))
Non-corrected
#using the clean culled df:
#average water year temperature
yearly_wy_aver_825 <- snotel_825_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp = mean(temperature_mean))
#Average temperature by day for all water years:
daily_wy_aver_825 <- yearly_wy_aver_825 %>%
group_by(daymonth) %>%
mutate(aver_day_temp = mean(aver_ann_temp))
#average mean temperature by day for the period of record:
daily_wy_aver_825 <- daily_wy_aver_825 %>%
group_by(daymonth) %>%
mutate(all_ave_temp = mean(daily_wy_aver_825$aver_day_temp))
# try to show all years as means.
daily_wy_aver2_825 <-daily_wy_aver_825 %>%
group_by(waterDay) %>%
mutate(date_temp = mean(temperature_mean))
daily_wy_aver2_825$date_temp <- signif(daily_wy_aver2_825$date_temp,3) #reduce the sig figs
ggplot(daily_wy_aver2_825, aes(x = waterDay, y = date_temp))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
standard_dev_825 <- daily_wy_aver_825 %>%
group_by(waterYear) %>%
mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_825 <- standard_dev_825 %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_825 <- standard_dev_all_825 %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_825 %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')
| waterYear | sd_2 |
|---|---|
| 1988 | 10.158111 |
| 1989 | 9.521283 |
| 1990 | 9.387663 |
| 1991 | 9.495098 |
| 1992 | 8.468342 |
| 1993 | 8.686860 |
| 1996 | 9.236157 |
| 2000 | 8.951568 |
| 2001 | 9.626977 |
| 2003 | 9.317887 |
| 2004 | 8.794074 |
| 2005 | 8.168991 |
| 2006 | 9.046952 |
| 2007 | 9.112464 |
| 2009 | 8.322586 |
| 2010 | 8.961897 |
| 2011 | 8.871929 |
| 2012 | 8.997194 |
| 2013 | 9.382743 |
| 2014 | 8.649080 |
| 2015 | 8.104451 |
| 2016 | 8.792359 |
| 2017 | 8.593127 |
| 2018 | 8.634239 |
| 2019 | 8.862025 |
| 2020 | 9.328094 |
| 2021 | 9.270956 |
| 2022 | 9.124264 |
ggplot(standard_dev_all_825, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Non-corrected standard deviation of SNOTEL 825 average temperatures for water years 2005-2021
sd_mk_825 <- mk.test(standard_dev_all_825$sd_2)
print(sd_mk_825)
##
## Mann-Kendall trend test
##
## data: standard_dev_all_825$sd_2
## z = -1.7188, n = 28, p-value = 0.08565
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## -88.0000000 2562.0000000 -0.2328042
sd_sens_825 <- sens.slope(standard_dev_all_825$sd_2)
print(sd_sens_825)
##
## Sen's slope
##
## data: standard_dev_all_825$sd_2
## z = -1.7188, n = 28, p-value = 0.08565
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.04047480 0.00485368
## sample estimates:
## Sen's slope
## -0.02125455
#using the clean culled df:
#average water year temperature
yearly_wy_aver_825_ad <- snotel_825_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_825_ad <- yearly_wy_aver_825_ad %>%
group_by(daymonth) %>%
mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_825_ad <- daily_wy_aver_825_ad %>%
group_by(daymonth) %>%
mutate(all_ave_temp_ad = mean(daily_wy_aver_825_ad$aver_day_temp_ad))
# try to show all years as means.
daily_wy_aver2_825_ad <-daily_wy_aver_825_ad %>%
group_by(waterDay) %>%
mutate(date_temp_ad = mean(temp_ad))
daily_wy_aver2_825_ad$date_temp_ad <- signif(daily_wy_aver2_825_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_825_ad, aes(x = waterDay, y = date_temp_ad))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
standard_dev_825_ad <- daily_wy_aver_825_ad %>%
group_by(waterYear) %>%
#filter(waterYear >= 1987 & waterYear <= 2021) %>%
mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_825_ad <- standard_dev_825_ad %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_825_ad <- standard_dev_all_825_ad %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_825_ad %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')
| waterYear | sd_2 |
|---|---|
| 1988 | 9.527196 |
| 1989 | 8.967107 |
| 1990 | 8.793798 |
| 1991 | 8.946793 |
| 1992 | 7.952054 |
| 1993 | 8.156645 |
| 1996 | 8.650702 |
| 2000 | 8.338506 |
| 2001 | 8.985410 |
| 2003 | 8.675745 |
| 2004 | 8.209448 |
| 2005 | 8.172186 |
| 2006 | 9.049660 |
| 2007 | 9.115791 |
| 2009 | 8.324666 |
| 2010 | 8.964371 |
| 2011 | 8.875058 |
| 2012 | 9.001462 |
| 2013 | 9.386059 |
| 2014 | 8.652436 |
| 2015 | 8.107386 |
| 2016 | 8.796286 |
| 2017 | 8.596142 |
| 2018 | 8.637340 |
| 2019 | 8.865259 |
| 2020 | 9.331765 |
| 2021 | 9.273785 |
| 2022 | 9.127176 |
ggplot(standard_dev_all_825_ad, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Morrisey-corrected standard deviation of SNOTEL 825 average temperatures for water years 1986-2021
sd_mk_825_ad <- mk.test(standard_dev_all_825_ad$sd_2)
print(sd_mk_825_ad)
##
## Mann-Kendall trend test
##
## data: standard_dev_all_825_ad$sd_2
## z = 1.0076, n = 28, p-value = 0.3137
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## 52.0000000 2562.0000000 0.1375661
sd_sens_825_ad <- sens.slope(standard_dev_all_825_ad$sd_2)
print(sd_sens_825_ad)
##
## Sen's slope
##
## data: standard_dev_all_825_ad$sd_2
## z = 1.0076, n = 28, p-value = 0.3137
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.009639429 0.034852064
## sample estimates:
## Sen's slope
## 0.01232876
Trapper Lake 827 Original 12/13/2004
Original 12/13/2004
snotel_827 <- SNOTEL_yampa_area %>%
filter(site_id == "827")
#str(snotel_827) # check the date, usually a character.
snotel_827$Date <- as.Date(snotel_827$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.
#THIS WILL CHANGE FOR EACH STATION
snotel_827_clean <- snotel_827 %>% # filter for the timeframe
filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
#filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers
addWaterYear() %>%
mutate(daymonth = format(as.Date(Date), "%d-%m")) %>%
na.omit()
#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))
snotel_827_clean <- snotel_827_clean %>%
group_by(waterYear)%>%
mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))
# Check for outliers
ggplot(snotel_827_clean, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
snotel_827_clean <- snotel_827_clean %>%
mutate(temp_diff = abs(temperature_min - temperature_max)) %>%
filter(temperature_mean > -40) %>%
filter(temperature_mean < 30)
ggplot(snotel_827_clean, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”
# filtering for temperature anomalies
#snotel_827_cull_count <- snotel_827_clean %>%
# filter(temperature_min > -40) %>%
# count(waterYear)
#snotel_827_cull_count
# filtering for too few observations in a year
snotel_827_cull_count_days <- snotel_827_clean %>%
group_by(waterYear) %>%
count(waterYear) %>%
filter(n < 350)
snotel_827_cull_count_days
## # A tibble: 18 x 2
## # Groups: waterYear [18]
## waterYear n
## <dbl> <int>
## 1 1986 331
## 2 1994 179
## 3 1995 323
## 4 1996 53
## 5 1997 278
## 6 1999 341
## 7 2001 348
## 8 2002 310
## 9 2003 340
## 10 2004 342
## 11 2009 342
## 12 2010 346
## 13 2011 326
## 14 2013 236
## 15 2014 241
## 16 2016 321
## 17 2017 219
## 18 2019 332
snotel_827_clean_culled <- snotel_827_clean %>%
filter(waterYear != "1986" & waterYear != "1994" & waterYear != "1995" & waterYear != "1996" & waterYear != "1997" & waterYear != "1999" & waterYear != "2001" & waterYear != "2002" & waterYear != "2003" & waterYear != "2004" & waterYear != "2009" & waterYear != "2010" & waterYear != "2011" & waterYear != "2013" & waterYear != "2014" & waterYear != "2016" & waterYear != "2017" & waterYear != "2019") #%>%
#filter(temperature_mean > -49)
ggplot(snotel_827_clean_culled, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
ggplot(snotel_827_clean_culled, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
temp_827_xts <- xts(snotel_827_clean_culled$temperature_mean, order.by = snotel_827_clean_culled$Date)
dygraph(temp_827_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)")
#snotel_827_clean_culled <- snotel_827_clean_culled %>%
# filter(temperature_mean > -30)
#temp_827_xts <- xts(snotel_827_clean_culled$temperature_mean, order.by = snotel_827_clean_culled$Date)
#dygraph(temp_827_xts) %>%
# dyAxis("y", label = "Daily mean temperature (°C)")
Original 12/13/2004
snotel_827_adjusted <- snotel_827_clean_culled %>%
mutate(temp_ad = if_else(Date < "2004-12-13", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))
Non-corrected
#using the clean culled df:
#average water year temperature
yearly_wy_aver_827 <- snotel_827_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp = mean(temperature_mean))
#Average temperature by day for all water years:
daily_wy_aver_827 <- yearly_wy_aver_827 %>%
group_by(daymonth) %>%
mutate(aver_day_temp = mean(aver_ann_temp))
#average mean temperature by day for the period of record:
daily_wy_aver_827 <- daily_wy_aver_827 %>%
group_by(daymonth) %>%
mutate(all_ave_temp = mean(daily_wy_aver_827$aver_day_temp))
# try to show all years as means.
daily_wy_aver2_827 <-daily_wy_aver_827 %>%
group_by(waterDay) %>%
mutate(date_temp = mean(temperature_mean))
daily_wy_aver2_827$date_temp <- signif(daily_wy_aver2_827$date_temp,3) #reduce the sig figs
ggplot(daily_wy_aver2_827, aes(x = waterDay, y = date_temp))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
standard_dev_827 <- daily_wy_aver_827 %>%
group_by(waterYear) %>%
mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_827 <- standard_dev_827 %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_827 <- standard_dev_all_827 %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_827 %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')
| waterYear | sd_2 |
|---|---|
| 1987 | 8.785352 |
| 1988 | 9.643636 |
| 1989 | 9.399426 |
| 1990 | 9.017517 |
| 1991 | 9.053817 |
| 1992 | 8.461780 |
| 1993 | 8.592203 |
| 1998 | 9.226580 |
| 2000 | 8.675610 |
| 2005 | 7.901837 |
| 2006 | 8.811208 |
| 2007 | 8.956871 |
| 2008 | 9.053611 |
| 2012 | 8.641372 |
| 2015 | 7.855779 |
| 2018 | 8.196461 |
| 2020 | 8.933495 |
| 2021 | 8.790312 |
| 2022 | 8.662496 |
ggplot(standard_dev_all_827, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Non-corrected standard deviation of SNOTEL 827 average temperatures for water years 2005-2021
sd_mk_827 <- mk.test(standard_dev_all_827$sd_2)
print(sd_mk_827)
##
## Mann-Kendall trend test
##
## data: standard_dev_all_827$sd_2
## z = -1.7493, n = 19, p-value = 0.08024
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## -51.0000000 817.0000000 -0.2982456
sd_sens_827 <- sens.slope(standard_dev_all_827$sd_2)
print(sd_sens_827)
##
## Sen's slope
##
## data: standard_dev_all_827$sd_2
## z = -1.7493, n = 19, p-value = 0.08024
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.06891402 0.00422489
## sample estimates:
## Sen's slope
## -0.03456914
#using the clean culled df:
#average water year temperature
yearly_wy_aver_827_ad <- snotel_827_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_827_ad <- yearly_wy_aver_827_ad %>%
group_by(daymonth) %>%
mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_827_ad <- daily_wy_aver_827_ad %>%
group_by(daymonth) %>%
mutate(all_ave_temp_ad = mean(daily_wy_aver_827_ad$aver_day_temp_ad))
# try to show all years as means.
daily_wy_aver2_827_ad <-daily_wy_aver_827_ad %>%
group_by(waterDay) %>%
mutate(date_temp_ad = mean(temp_ad))
daily_wy_aver2_827_ad$date_temp_ad <- signif(daily_wy_aver2_827_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_827_ad, aes(x = waterDay, y = date_temp_ad))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
standard_dev_827_ad <- daily_wy_aver_827_ad %>%
group_by(waterYear) %>%
#filter(waterYear >= 1987 & waterYear <= 2021) %>%
mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_827_ad <- standard_dev_827_ad %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_827_ad <- standard_dev_all_827_ad %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_827_ad %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')
| waterYear | sd_2 |
|---|---|
| 1987 | 8.205336 |
| 1988 | 9.015708 |
| 1989 | 8.826659 |
| 1990 | 8.399945 |
| 1991 | 8.491804 |
| 1992 | 7.896678 |
| 1993 | 8.029814 |
| 1998 | 8.583994 |
| 2000 | 8.036502 |
| 2005 | 7.717045 |
| 2006 | 8.811366 |
| 2007 | 8.957047 |
| 2008 | 9.053812 |
| 2012 | 8.641373 |
| 2015 | 7.856092 |
| 2018 | 8.196931 |
| 2020 | 8.933512 |
| 2021 | 8.790758 |
| 2022 | 8.663052 |
ggplot(standard_dev_all_827_ad, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Morrisey-corrected standard deviation of SNOTEL 827 average temperatures for water years 1986-2021
sd_mk_827_ad <- mk.test(standard_dev_all_827_ad$sd_2)
print(sd_mk_827_ad)
##
## Mann-Kendall trend test
##
## data: standard_dev_all_827_ad$sd_2
## z = 0.41983, n = 19, p-value = 0.6746
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## 13.00000000 817.00000000 0.07602339
sd_sens_827_ad <- sens.slope(standard_dev_all_827_ad$sd_2)
print(sd_sens_827_ad)
##
## Sen's slope
##
## data: standard_dev_all_827_ad$sd_2
## z = 0.41983, n = 19, p-value = 0.6746
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.03007505 0.05895184
## sample estimates:
## Sen's slope
## 0.012232
Note: figure captions are incorrect.